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ABSTRACT 

In this study we demonstrate for the first time that the unified Monte Carlo 
approach can be applied to model gas-grain chemistry in large reaction networks. 
Specifically, we build a time- dependent gas-grain chemical model of the interstel- 
lar medium, involving about 6000 gas-phase and 200 grain surface reactions. 
This model is used to test the validity of the standard and modified rate equa- 
tion methods in models of dense and translucent molecular clouds and to specify 
under which conditions the use of the stochastic approach is desirable. 

Two cases are considered: (1) the surface mobility of all species is due to 
thermal hopping, (2) in addition to thermal hopping, temperature-independent 
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quantum tunneling for H and H2 is allowed. The physical conditions characteris- 
tic for the core and the outer region of TMCl cloud are adopted. The gas-phase 
rate file RATE 06 together with an extended set of gas-grain and surface reactions 
is utilized. 

We found that at temperatures 25-30 K gas-phase abundances of H2O, NH3, 
CO and many other gas-phase and surface species in the stochastic model differ 
from those in the deterministic models by more than an order of magnitude, at 
least, when tunneling is accounted for and/or diffusion energies are 3x lower than 
the binding energies. In this case, surface reactions, involving light species, pro- 
ceed faster than accretion of the same species. In contrast, in the model without 
tunneling and with high binding energies, when the typical timescale of a surface 
recombination is greater than the timescale of accretion onto the grain, we ob- 
tain almost perfect agreement between results of Monte Carlo and deterministic 
calculations in the same temperature range. At lower temperatures (~ 10 K) 
gaseous and, in particular, surface abundances of most important molecules are 
not much affected by stochastic processes. 

Subject headings: astrochemistry ISM: abundances, clouds, molecules — molec- 
ular processes - radio lines: ISM — stars: formation 



1. Introduction 

Chemical processes are discrete in nature, and it has been realized long time ago that 
their microphy sically correct theoretical treatment should rest upon stochastic methods 



( lGillespielll976l ). While in cellular biology Monte Carlo simulations of chemical processes in 



cells are widespread, most models in astrochemistry are based on deterministic rate equa- 



"cell-like" obiects like dust grains and PAHs 


(Tielens & Haeen 


1982; 


Caselli et al. 


1998; 


Herbst & Shematovich 


2003; 


Lipshtat & Biham 


2003; 


Stantcheva & Herbst 


2004 


). The con- 



cept of mass-action kinetics, which hes behind rate equations, being adequate for gas-phase 
reactions, is not appropriate in the case of small surface population (number of particles 
< 1 per grain). This situation is often the case in grain surface catalysis. Rates of surface 
reactions obtained through this approach may be strongly incorrect, leading to an improper 
estimation of timescales of fundamental processes, like H2 formation, and significant errors 
in abundances of other species. 



Up to now, a number of attempts have been made to assess the significance of stochastic 
effects for astrochemical simulations and to develop suitable numerical methods. In the 
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context of ast rochemistry first a ttemp t to treat grain surface cliemistry stocliastically lias 
been made by iTielens fc HagenI (jl982l ). In tliis study, equilibrium abundances of gas-phase 
species were used to obtain accretion rates and to calculate time-dependent populations of 
surface species through Monte Carlo approach. All surface species but H2 were assumed to 
remain on the grain surfaces. 

Subsequent attempts to account for the discrete nature of surface reactions can be di- 
vided into two categories. The first ca tegory encornpasses different modifications to th e rate 
equations for grain surface chemistry (ICaselli et al.lll998l . l2002l : IStantcheva et al.ll200ll ). The 
main idea behind these modifications is to restrict rates of diffusive surface reactions by ac- 
cretion or desorption rates of reactants. Thus, the so-called accretion-limited regime is taken 
into account, when the rate of a particular surface reaction is determined by the flux of reac- 
tants accreting on a grain and not by their diffusion over the surface. For simple reacting sys- 
tems, results obtained with this semi-empirical approach were found, at least, to be in closer 
agreement with Monte Carlo (MC) computations than the treatment by rate equations alone. 
However, the applicability of the modified rate equations (MRE) to large chemical networks 
has not been tested so far. Moreover, it is impossible to investigate purely stochastic effects. 



like bist ability, with this technique , still purely deterministic in its nature (ILe Bourlot et al. 



I993I : IShalabiea fc GreenberelEggsl : llee et al.lllQQsl : ISoger fc SternberelboOGf ) 



In the second category of studies the master equation is solved with different meth- 
ods. Because the chemical master equation describes probabilities of a chemically reacting 
system to be in all possible states of its phase space and the number of these states grows 
exponentially with the number of species, their direct integration is only possible for ex - 
tremely simple syst ems, consisting of a v ery fe w different components. iBiham et al.l (120011 ). 
Green et al.l (I2OOII ) and iLipshtat et al.l (120041 ) have performed a direct integration of the 
master equation to investigate H2 formation. These studies demonstrated that in some cases 
th is fundamental pr o cess c a nnot be adequately desc ribed with the rate equatiori approach. 
In Istantcheva et al.l J2002h . ICaselh et all J2002h and Istantcheva fc HerbstI J2003h direct in- 
tegration of the master equation was applied to study the evolution of a simple H — O — CO 
surface chemical network and deuterium fractionation. Direct integration of the master 
equation in more complex networks of grain surface reactions is hampered by the lack of 
appropriate computing power. To the best of our knowledge, the only successful simulation 
of a more extended surface chemical network in astrochemistry, consi sting of 19 reactions. 



with direct integration of the master equation was accomphshed by IStantcheva fc Herbst 



A Monte Carlo approac h based on continuous-time random-walk method has been devel- 
oped by lChang et al.l (|2005l ) to model recombination of hydrogen atoms on interstellar grains. 
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The advantage of this method is that it allows to model formation of H2 a nd other molecules 
on sur faces of arbitrary roughness, which has later been demonstrated by ICuppen fc Herbst 
J2007h . 



To relax computational requirements n eeded to solve the raaster equat ion, the moment 
equati on approximation was suggested by iLipshtat fc Bihanj (120031 ) and iBarzel fc Biham 
(120071 ) ■ Currently, this is the most promising approach to efficient simulations of gas-grain 
systems with complex grain surface chemistry, because the moment equations can be easily 
combined with rate equations used to simulate gas-phase chemical processes. But, strictly 
speaking, the moment equations are an approximation and their validity should be verified 
by a comparison with exact methods. 

The only feasible technique allowing us to obtain an exact solution to the mas ter equa- 
tion fo r complex chemically reacting systems is Monte Carlo method developed by I Gillespie 
(119761 ). This approach is widely used in molecular biology for simulation of chemical pro- 
cesses in living cells and can also be applied to astrochemical problems. Gillespie's Stochastic 
Simul ation Algorithm (SSA) was first used in complex astrochemical networks by ICharnley 
(119981 ) . where SSA was applied to gas-phase chemistry onl y. Later a simple g rain surface 
H — O — CO chemical mo del, similar to the one considered in lCaselli et al.l (119981 ) . was devel- 
oped by lCharnleyl (I2OOII ). In this model, gas-phase abundances of species were assumed to be 
constant and desorption was neglected. No coupled gas-grain chemical model was studied. 

The chemical master equation is usually considered as a method to model surface reac- 
tions. However, due to its universality, it can be used to model the gas-phase chemistry as 
well. The reason why the rate equations are so popular as a tool to study gas-phase chem- 
istry in astronomical environments is the fact that they are easier to implement and much 
more computationally effective. However, if one uses the rate equations to model gas-phase 
chemistry and some stochastic approach to model surface chemistry one has to find a way 
to match the two very different computational techniques. The first study in which a Monte 
Carlo t reatment of g r ain su rface chemistry was coupled to time-dependent gas-phase chem- 



istry is IChang et al.l (120071 ). In this paper iterative technique is utilized to combine very 
different Monte Carlo and rate equations approaches. Calculations are performed for the 
time interval of 2 ■ 10^ years and a relatively simple H — O — CO surface network consisting 
of 12 reactions is used. 

Even though remarkable progress has been made to develop stochastic approaches to 
astrochemistry, up to now there are no studies in which a stochastic treatment of a complex 
grain surface chemical network is fully coupled to time-dependent gas-phase chemistry. The 
situation gets less complicated when the same method is used for both gas-phase and surface 
chemistry. In this paper, for the first time, we present a "complete" time-dependent chemical 
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gas-grain model, calculated with Monte Carlo approach. We solve the master equation 
with a Monte Carlo technique, using it as a single method to treat gas-phase reactions, 
surface reactions, and gas-grain exchange processes simultaneously. The gas-phase chemical 
network includes more than 6000 reactions while surface network consists of more than 200 
reactions. This Monte Carlo model is used to test the validity of rate equations and modified 
rate equations for a range of physical conditions, typical of diffuse and dense clouds. So far, 
no comparison of RE and MRE techniques against rigorous stochastic approach over a wide 
range of physical conditions, typical for astrophysical objects, has been made. 

There is clear observational evidence that interstellar grains have complex mineral com- 
position and may have either amorphous or crystalline structure. In general, it is believed 
that interstellar dust parti cles are made of oliv ine-like silicates and some form of amorphous 
carbon or even graphite (jPraine &: Led Il984l ). The location on a surface which may be 
occupied by an accreted molecule is commonly referred to as a surface site. The density 
of sites and binding energies for specific species strongly depend on the dust material and 
its structure. Bi nding energies E y, ^ that is, potential barriers between adjacent sites, have 
been assumed bv iHasegawa et a.1.1 f 1992[) to be roughly 0.3 o f corr espo nding desorpt i on en- 



ergies Eq. In IPirronello et all (ll997al Jbl. Il999[ ). iKatz et all (Il999f ) and IPerets et all (120051 ) 



experimental results on the H2 formation, the ratio of diffusion energy to desorption energy, 
i^bZ-^D, was found to be close to 0.77 for H atoms. This ratio determines the surface mobility 
of species and is of utter importance for estimating surface reaction rates. 

Another factor defining rates of reactions involving H and H2 (and their isotopologues) is 
related to the possibility of their tunneling. Non-thermal quantum tunneling of these species 
allows them to scan grain surface and find a reacting partner quickly. This mechanism has 
been suggested as an explanation o f efficient forrn ation in the ISM. Later analysis of the 
cited experimental results has lead iKatz et al.l (119991 ) to the conclusion that tunneling does 
not happen. The absence of tunneling implies lower mobility, so that surface processes are 
essentially reduced to accretion and deso rption. However, this co nclusion turned out to be 



not the final one. In other studies (e.g. ICazaux fc Tielenj |200J) alternative explanations 



of these expe riments hav e been proposed and some shortcomings of the theoretical analysis 
performed in iKatz et al.l (119991 ) have been found. So, up to now the question of presence 
or absence of tunneling effects on grain surfaces is not settled, and the E\,/Ej) ratio is not 
well constrained. With this in mind, two different models of grain surface reactions are 
considered: (1) the surface mobility of all species is caused by thermal hopping only with 
high El,/ Ed ratio equal to 0.77, (2) temperature-independent quantum tunneling for H and 
H2 is included in addition to thermal hopping with low E\^/Eiy ratio equal to 0.3. 



The organization of the paper is as follows. In Section 2 we describe physical conditions 
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adopted in simulations and chemical model used. Section 3 contains basics of stochastic 
reaction kinetics and detailed description of the Monte Carlo code used for simulations. 
In Section 4 the validity of rate equations and modified rate equations is checked against 
Monte Carlo method. First, global agreement between methods is investigated. Then, some 
interesting species are discussed separately. In Section 5 a discussion is presented. Section 6 
contains the conclusions. 



2. Modeling 
2.1. Physical conditions 



In the present study, we consider grain surface reactions under physical conditions typ- 
ical of irradiated translucent clouds, cold dark cores and infrared dark clouds. We consider 
temperatures T between 10 and 50 K, densities n(H) between 200 and 2 ■ 10^ cm~^ (three 
values of n(H) and five values of T), and visual extinctions of 0.2 mag at n CH) = 200 cm~ ^, 



2.0 rnag at njB) = 2 • 1 0^ cm ^ and 15 mag at n(H) = 2 ■ 10^ cm ^ (e.g., ISnow fc McCall 



20061 : iHassel et al.ll2008l ). This roughly corresponds to the distance of order of 0.4-0.5 pc 
from the cloud boundary (at constant density). The cloud is illuminated by the mean inter- 
stellar diffuse UV field. The dust temperature is assumed to be equal to gas temperature. 
We do not study the earliest stage of the molecular cloud evolution, which is essentially the 
stage when atomic hydrogen is converted into hydrogen molecules. The chemical evolution 
is simulated assuming static physical conditions. While it is still a matter of debate how 
long does it take for a typical isolated cloud core to become gravitationally unstable and to 
start collapsing, we as sume a short timescale of 1 Myr which see ms to be appropriate for the 
TMCl cloud (see e.g. [Roberts et aPbooi ISemenov et"aDl2004[ ). 



2.2. Chemical model 



The chemical model is the same as described in IVasyunin et al.l (120081 ). The gas-phase 
reactions and their rates are taken from the RATE 6 database, in which the effects of dipole- 
enhanced ion- neutral rates are taken into account (jWoodall et al.l 120071 ). All reactions with 
large negative activation barriers are excluded. The rates of photodissoci ation and photoion 
i zation of molecular species by interstellar UV photons are taken from Ivan Dishoeck et al. 
(120061 ). The self- and mutual shielding of CO ari d H2 against UV photodissociation are com- 
puted as des cribed in Ivan Zadelhoff et al.l (120031 ) using pre-calculated factors from Tables 10 
and 11 from iLee et al.l (119961 ). Ionization and dissociation by cosmic ray particles are also 
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considered, with a cosmic-ray ionization rate of 1.3 -10 s ^ (ISpitzer fc Tomaskolll968l ) 



The gas-grain interactions include accretion of neutral species onto dust grains, thermal 
and photodesorption of mantle materials, dissociative recombination of ions on charged grains 
and grain re-charging processes. The dust grai ns are uniform 0.1 nm spherical particles made 
of amorphous silicates of olivine stoichiometry (jSemenov et al. with a dust-to-gas mass 

ratio of 1%. The s ticking probabihty is 100 %. Desorption energies and surface reaction 
list are taken from iGarrod fc HerbstI (j2006( ). 



The grain surface is assumed t o be compact, with surface density of 2 ■ 10^^ sites cm ^, 



which gives ~ 3 ■ 10^ sites per grain (IBiham et al.ll200ll ). We employ two models to calculate 
the rates of surface reactions (Tabled]). In Model T (T for tunneling) the tunneling timescale 
for a light atom to overcome t he potential barr i er ari d migrate to another potential well is 
computed using Eq. (10) from iHasegawa et al.l ( 19921 ). with the barrier thickness of lA. In 
Model H (H for hopping) we do not allow H and H2 to scan surface sites by tunneling. 
The diffusion timescale for a molecule is calculated as the timescale of thermal hopping 
multiplied by the tota l number of surface site and is given by Eq. (2) and Eq. (4) from 
Hasegawa et al.l (119921 ). The hopping rates are sensitive to the adopted values of diffusion 
energy, which are not well constrained. Thus, we consider low and high diffusion energies 
which are calculated from ad opted desorption ene rgies by multiplyi ng them by factors of 0.3 
(like in IHasegawa et al.lll992h (Model T) and 0.77 JXatz et al.lll999h (Model H), respectively. 
For all considered models, the total rate of a surface reaction is calculated as a sum of diffusion 
or tunneling rates divided by the grain number density and multiplied by the probability of 
reactions (100% for processes without activation energy). 

Overall, our network consists of 422 gas-phase and 157 surface species made of 13 
elements, and 6002 reactions inclu ding 216 s urface reactions. As initial abundances, we 
utilize the "low metallicity" set of iLee et al.l (119981 ). where abundances of heavy elements 
in the gas are assumed to be severely depleted. All hydrogen is molecular initially. The 
chemical evolution for 1 Myr in the classical deterministic approaches is computed with the 
fa st "ALCHEMIC" co djll in which the modified rate equations are implemented according 
to Caselli et al.l (1998). No fur ther modification for reactions with activation energy barriers 



is used (see lCaselli et al.ll2002l ) 



^Available upon request: semenov@nipia.de 
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3. Stochastic reaction kinetics 
3.1. Theoretical foundations 

We consider a chemically reacting system which consists of N d ifferent types of species 



{5*1 ... Sn} and M chemical reactions {Ri ... Rm} (iGillespid 119761 ). All these species are 
contained in a constant volume Q, in which local thermal equilibrium is reached, so that the 
system is well mixed. We denote the number of species of type i at time t as Xi(t). The 
ultimate goal is to determine the state vector X{t) = {Xi{t)...Xfyf{t)} of the system at any 
given time t > to, assuming certain initial conditions X(to) = -^o- 

Let us assume that each reaction Rj retains properties of a Markov chain and can be 
considered as a set of independent instantaneous events. Each chemical reaction is described 
by two quantities, the discrete state vector z7 = {uij, unj} and the propensity function aj. 
The components of the state vector represent the net change in populations of the ith 
species due to the jth reaction. By definition, the propensity function 

a,iX)dt (1) 

is the probability for a reaction Rj to occur in the volume Q over the time interval dt. 
In analogy to the reaction rate in the rate equation approach, for a one-body reaction the 
propensity function aj can be expressed through rate constants and numbers of reactants as 

aj{X) = CjXi. (2) 

Here, Cj is the rate constant of the jth reaction and Xi is the absolute population of reagent. 
For a two-body heterogeneous reaction this expression changes as follows: 

aj{X) = CjXiX2, (3) 

where cj is the reaction rate constant, xi and X2 are the abundances of first and second 
reactants. If we deal with a two-body homogeneous reaction (like H -|- H — > II2), the 
expression for its propensity function becomes somewhat more complicated: 

a,(X) = |xi(xi-l) (4) 

This expression refiects the fact that the rate of a homogeneous two-body reaction is propor- 
tional to the number of all possible pairs of its reactants. The term [xi — 1) is important for a 
proper calculation of the reaction rate in the stochastic regime when the average reagent pop- 
ulation is close to unity and cannot be properly reproduced by the rate equation approach. 
Note that the abundances of species are always integer numbers in these expressions. 
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These definitions allow to characterize the microphysical nature of chemical processes 
and to establish a basis for stochastic chemical kinetics. A vector equation that describes 
the temporal evolution of a chemically reacting system by stochastic chemical kinetics is the 
chemical master equation. To derive this equation one has to introduce the probability of a 
system to be in the state X(t) at a time t > to: 



P{X,t I Xo,to) = Proh{X{t) = X I X{to) = Xo 



(5) 



Here P{X,t \ Xo,to) is the conditional probability density function of the time-dependent 
value X. Temporal changes of X are caused by chemical reactions with rates defined by 
Eq. ([T]). Therefore, the change of P{X,t \ Xo,to) over the time interval dt is the sum of 
probabilities of all possible transitions from the state Xq to the state X{t): 

P(x,t + dt\ Xo,^o 



P \ X,t I Xq, to 



M 

1 - J2 {cij{^)dt^ 



+ 



M 

- zT;-, t I Xo, to) ■ (ai(X - Uj)dt^ 



j=i 

The first term in this equation is derived from the fact that the system is already in the state 
X and describes the probability for the system to leave the state during the time interval 
dt. The second term defines the total probability for the system to reach the state X from 
states X — Vj during the time interval dt due to a reaction Rj. Finally, we formulate the 
chemical master equation: 

M 

^K(X - u,)P{X - i7^,t I Xo,to) - a,(X)P(X,t | Xo,to)] (6) 



dP{X,t\Xo,to) 



dt 



This vector equation looks similar to the conventional set of balance equations. However, it 
properly takes into account the discrete nature of chemical processes and, thus, has a much 
wider range of applicability. Unfortunately, an analytic integration of the chemical master 
equatio n is only possible in a very limited number of cases (several species and reactions. 



see e.g. 
to solve this equation. 



Stantcheva fc HerbstI |2004| ) . Thus, one has to rely on various numerical techniques 



3.2. Implementation of the Monte Carlo algorithm 



In principle, any normalization of abundances can be used in Eq. ([6]). In particular, 
if we take to be a unit volume, we end up with the usual number densities, which are 



- 10 - 



widely used in astrochemical models. When dealing with mixed gas-dust chemistry, one has 
to be more careful with the normalization. As surface reactions can only proceed when both 
reactants reside on a surface of the same grain, we define Q as the volume of the interstellar 
medium that contains exactly one dust grain. 



"MC — 5 [I ) 

Pgas7 



where r is the grain radius, pdust is the mass density of grain material, pgas is the mass gas 
density, and 7 is the dust-to-gas (0.01) mass ratio. Here grains are assumed to be spheres of 
equal size. As we also assume that the interstellar medium is well-mixed, so that the volume 
Que is representative for any volume of the real interstellar medium with the same physical 
conditions. 

This setup allows to construct the unified master equation for all three kinds of processes 
listed above (gas phase reactions, accretion/desorption processes and grain surface r eactions) . 



This e quation is solved with the stochastic Monte Carlo algorithm described in I Gillespie 



( 1l976l ). We implemented this technique in a FORTRAN77 code which allows to simulate 
all chemical reactions in the network in a self-consistent manner. Such a "brute-force" 
approach requires substantial CPU power and cannot be utilized in massive calculations. 
In the present study, it is used as a benchmark method to simulate chemical evolution in 
astrophysical objects. A typical run on a single Xeon 3.0GHz CPU takes between 10 hours 
and several days of computational time, and involves several billions of time steps. If one 
is only interested in a single point then the model can be used not only for benchmarking, 
but also for practical purposes. Unfortunately, it becomes impractical if abundances are 
computed for a number of spatial locations. High density also slows down the computation 
significantly. 

Due to the fact than the Monte Carlo technique operates with integers, the smallest 
abundance of a molecule is 1. Given the dust properties in our model (radius r = 10~^ cm, 
density pdust = 3 g cm~^), a dust-to-gas ratio of 1/100 means that volume per one grain is 
about 10^^/n(H) cm^, where n(H) is the number density of hydrogen nuclei. Therefore, an 
absolute population of 1 corresponds to a relative abundance of 10~^^ with respect to the 
total number of hydrogen atoms (see Eq. [7]). This is the lowest abundance directly resolvable 
in our model. However, the huge amount of tiny time steps taken during the calculations 
allows to average stochastic abundances over wider time spans of 10-1 000 years and push 
the smallest abundance resolvable by the Monte Carlo code below 10^^^: 
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where Xi is the abundance of species after the ith time step, the Ati is the time step in 
seconds and is the amount of time steps over which the averaging is performed. In such 
an averaging procedure, the noise in the random variable decreases with N as and 
enables us to resolve abundances below 10~^^. In our simulations we averaged abundances 
over 10^ time steps. Therefore, the resolution of our calculations is about 10~^^ — 10~^^ 
cm~^. The average value of the Monte Carlo time step in Gillespie's algorithm is equal to 
an inverse sum of rates of all chemical reactions in the model. For example, in the dense 
cloud model many molecules form rapidly at initial times and less actively toward the end of 
evolution (Fig. [1]). Time intervals, corresponding to 10^ time steps in a medium with density 
10^ cm~^ vary from 10 to 7 ■ 10^ years for early and late time moments, correspondingly. 



4. Comparison of the Monte Carlo method with rate equations 



While the Monte Carlo (MC) technique seems to be the most adequate method to solve 
the chemical master equation, it is computationally extremely demanding. Rate equations 
(RE) do not account for the stochastic nature of surface processes and thus may produce 
spurious results in some circumstances. On the other hand, they do require much less 
computer power and are easier to handle. Therefore rate equations will remain an important 
tool in theoretical astrochemistry. 

Given the importance of surface processes, it is thus necessary to isolate regions in the pa- 
rameter space where rate equations should not be used at all or, at least, should be used with 
caution. Similar comparisons , having been made so far, are based on a very limited number 



of species and reactions (e.g. ICharn 



Case! 



20031: 



i et 



aDl200^; 



Stantcheva et al. 



Barzel fc BihamlboOTl : IChang et al.l 



1998 



20021: 



2001 



Ruffle fc Herbst"2000l: ICreen et al. 



Stantcheva fc Herbst 



2003i.,2004 



2001 



Lipshtat fc Biham 



200 



3). As we implemented a unified Monte Carlo 



approach to the gas-phase and surface chemistry, which is capable to treat a "regular size" 
chemical network, we are able to perform a comprehensive comparison between predictions 
of the MC and the (M)RE approaches. 

For both the Model T and the Model H, the rate equations (RE), the modified rate 
approach (MRE) (for surface chemistry) and the Monte Carlo method are used to simulate 
gas-phase and surface chemistry. The differences in abundances discussed are related to 
surface chemistry and/or dust-gas interactions. For the gas-phase chemistry, all models 
produce identical results. Each run represents a combination of the tunneling treatment (T 
or H) and the surface chemistry treatment (RE, MRE, or MC). In total, we considered six 
models denoted as T-RE, T-MRE, etc. 
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4.1. Global Agreement 

To perform an overall comparison, in Figures [2] and [3] we show diagrams of global 
agreement between MC, RE, and MRE calculations as a function of n and T for evolutionary 
times 10^, 10^, and 10® years. For a given species, the two methods, being compared, are 
assumed to agree at time t if they produce abundances, that differ by no more than an 
order of magnitude. As we have mentioned above, in the MC method we are able to resolve 
abundances down to 10~^^, which would lead to a formal disagreement between abundances, 
say, 10~^^ in the MC method and 10"^® (or smaller) in the RE method (all abundances are 
relative to the total number of H atoms). This kind of disagreement is not meaningful from 
the observational point of view because such small relative abundances are extremely difficult 
to observe with a satisfactory accuracy. In addition, even though lower abundances can be 
calculated by the RE method, their actual accuracy is limited by the numerical interpolation. 
In the following discussion, species, for which both methods predict relative abundances less 
than 10^^^, are excluded from further consideration. In Figures [2] and [31 contours labeled 
with 0.9 mean that the two methods give abundances that differ by less than an order of 
magnitude for 90% of species with abundances higher than 10~^^ in both methods. 

An order-of-magnitude agreement criterion may seem to be too coarse. However, when 
stochastic effects do not play an important role, all three method produce results, which 
are nearly identical. To illustrate this, we show in Figure H] the evolution of the H2CCCC 
abundance at n = 2 ■ 10^ cm~^ and T = 10 K, simulated with the MC and RE methods. The 
species is chosen because its abundance strongly varies with time. It can be seen that the 
difference between the two methods is very small. 

4.1.1. Model T 

In Model T the species are highly mobile — atomic and molecular hydrogen because of 
tunneling, and all the other species because of the low energy barrier for diffusion. Mobility 
drives rapid surface reactions with rates which sometimes exceed the accretion/desorption 
rate. Surface reactions, which are on average faster than accretion, are not permitted in 
the MC approach when surface populations of their reagents are low (~ 1), but they are 
allowed in the RE approach, where the surface abundance of a species can be much less than 
one. Because of that, differences between Monte Carlo and rate equation methods are quite 
noticeable in Model T. 

At temperatures below 20 K and above ~ 40K the agreement between the stochastic 
method and the rate equation method (both RE and MRE) is about 85% or better at 
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all times and densities (Figure [2]). Around T ~ lOK the residual discrepancy is mainly 
caused by complex surface species, which have zero abundances in Model T-MC and are 
overproduced in Model T-RE and Model T-MRE, so that their abundances are just above 
the adopted cutoff of 10~^^. If we would raise the cutoff to 10~^^, the agreement would be 
almost 100%. At T ~ 50K the agreement is also nearly perfect, as the chemistry is almost 
a purely gas-phase chemistry under these conditions (dust temperature being equal to gas 
temperature). 

All the major discrepancies are concentrated in the temperature range between 25- 
30 K. In this range the accretion rate (which depends on gas temperature) is high enough to 
allow accumulation and some processing of surface species, while the correspondingly high 
desorption rate (which depends on dust temperature and Ay) precludes surface production 
of complex molecules. The latter process is adequately described in the MC runs only. 

At earlier times {t < 10^ years) the largest differences in this temperature range are 
observed at lower densities. They are caused by an overproduction of 'terminal' surface 
species, like water, ammonia, hydrogen peroxide, carbon dioxide, etc., in the RE and MRE 
calculations. At later times (t ~ 10^ years) situation changes drastically. The discrepancies 
shift toward higher density, and the overall agreement falls off to about 60% for both RE 
and MRE calculations. Even though the stochastic chemistry is only important on dust 
surfaces, by this time it also influences many gas-phase abundances due to effective gas- 
grain interactions. 

At if: = 10^ years, in the RE calculation at T = 30 K, among 234 species with abundances 
higher than 10~^^, 95 species (including 81 gas-phase species) have an order of magnitude 
disagreement with the MC run. In the MRE run at T = 25 K, the number of species with 
abundances above 10~^^ is 221. Among them 83 species disagree with the MC results by 
more than an order of magnitude. In both calculations these are mostly carbon-bearing 
species, like CO, CO2, CS, to name a few, and, in particular, carbon chains (like some 
cyanopolyynes) . It is noteworthy that not only the abundances of neutral species disagree, 
but also the abundances of some ions, including C^ and 

4.1.2. Model H 

In Model H the mobility of all species on the surface (including H and H2) is caused 
by thermal hopping only, which is slow because of the high ^^bZ-E-D ratio. One of the con- 
sequences of these slow surface reaction rates is that results produced by the modified rate 
equations almost do not differ from results of conventional rate equations. The essence of 
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the modification is to artificially slow down surface reactions to make their rates consistent 
with accretion/desorption rates. In Model H all surface reactions are slow anyway, so that 
modifications never occur. Thus, we only compare the results of the MC and RE calculations. 

At t = 10^ years the percentage agreement between the two runs never falls below 
80% and most often is actually much better than the adopted order-of-magnitude criterion. 
At later times, the discrepancy appears in the same temperature range as in Model T, 
i.e., between 25 and 30 K. Unlike Model T, both at t = 10^ years and t = 10^ years the 
disagreement shows up at low density. The set of discrepant species is nearly identical at 
both times and consist mostly of ices, overproduced by rate equations. These ices include 
some key observed species. For example, NH3, H2O, CO2 ices are all abundant in the 
Model H-RE despite the low density, with a surface H2O abundance reaching 3 ■ 10^'' by 10^ 
years. On the other hand, the Model H-MC results in very small or zero abundances for 
the same species. More chemically rich ices in Model H-RE bring about enhanced gas-phase 
abundances for the same molecules. The only species which is underproduced by the rate 
equation approach is surface C2O. 

4.2. Selected Species 

The diagrams presented above indicate areas in the parameter space, relevant to molec- 
ular clouds, where deterministic methods fail to describe stochastic surface processes. How- 
ever, a global agreement does not necessarily imply that the abundances of key species are 
also correctly calculated. In the following we discuss a few important gas-phase species and 
ices, for which an order of magnitude (or more) disagreement has been found. 

4.2.1. Model T 

Because surface reaction rates are high in this model, there are many species for which 
stochastic and deterministic methods give quite different results. The most abundant trace 
molecule, CO, is among these species. In Figure we show agreement diagrams for Mod- 
els T-RE and T-MRE (top row) and the CO abundance evolution for n — T combinations 
corresponding to the worst agreement. Both the RE and MRE methods fail to reproduce 
the late CO evolution at high density and temperatures of 25-30 K (where the global agree- 
ment is also worst in this model), but each in a different way. While in Model T-RE CO 
is underabundant by an order of magnitude, in Model T-MRE it is overabundant by the 
similar amount. This difference is related to the treatment of CO — > CO2 conversion on 
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dust surfaces, where underabundance of the CO ice in Model T-RE and overabundance in 
Model T-MRE are observed at almost all times. At later time, these surface abundances 
just start to propagate to gas-phase abundances. 

Overproduction of the CO2 ice in Model T-RE also consumes carbon atoms which would 
otherwise be available not only for CO, but also for more complex molecules, in particular, 
carbon chains, starting from C2, observed in diffuse clouds, and C2S, used as a diagnostics 
in prestellar cores. In the RE and MRE calculations these molecules, like CO, show trends 
opposite to that of CO2 ice: when surface CO2 is overabundant with respect to the MC 
model, carbon chains are underabundant and vice versa. 

When CO molecule sticks to a grain, it either desorbs back to the gas-phase, where it 
may participate in further processing, or is converted into CO2. As the desorption energy 
is quite high for CO2, it acts as a sink for carbon atoms at these temperatures. In the 
accretion-limited regime, rate equations overproduce CO2 with respect to the MC method. 
Thus, we see less CO and other carbon-bearing species in the gas-phase. The MRE method 
helps, probably, too much to account for accretion-limited CO processing and essentially 
quenches surface CO — > CO2 conversion. It is worth noting that at high density in this 
particular temperature range (and for adopted CO desorption energy) CO balances between 
complete freeze-out (at T < 20K) and near absence on dust surface (at T > 30 K). Thus, 
even relatively minor changes in treatment may lead to noticeable consequences. 

Another carbon-bearing molecule that shows an interesting behavior is methanol in 
the gas phase. Diagrams, comparing its abundance computed with the MC approach and 
with the rate equations, are presented in Figure [61 The agreement between results of the 
MC calculations and the RE calculations is almost perfect, at least, within the scatter 
produced by the MC simulations. However, the modified rate equations, which tend to 
improve agreement between stochastic and deterministic methods for many other species, in 
this particular case underpredict methanol abundance by more than two orders of magnitude 
in comparison with the MC run. The same is true for surface methanol as well. It is 
interesting to note that the region in the parameter space, where MC and MRE methanol 
abundances disagree, extends down to lowest temperatures considered. 

An interesting pattern is presented by the abundance of water, ammonia and their ices 
(Figure [7]). These species are severely overproduced by the RE methods at t = 10^ years, 
but after this time the modified rate equations are able to restore the agreement with the 
MC method, so that by 10^ years difference between MC and MRE runs is not significant. 
This is an example of a situation where the MRE help to produce realistic results, at least, 
at later times. We show plots only for gas-phase abundances because these are what is really 
observed. However, these abundances are only a reflection of surface processes, and surface 
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water and ammonia abundances behave in exactly the same way. It is hard to name a single 
reaction which is responsible for the difference between results of RE and MRE methods in 
this case. Let's consider water as its chemistry is somewhat simpler. The primary formation 
reaction for water is H + OH, with only a minor contribution from H2O2 (see below). 
However, hydroxyl is not produced in H + O reaction, as one might have expected. Because 
of the paucity of H atoms on the grain surface, hydroxyl formation is dominated by reaction 

O + HNO — ^ NO + OH. 

Abundance of HNO is restored in reaction 

H + NO — ^HNO. 

These two reactions form a semi-closed loop for which the sole result is OH synthesis out 
of an H atom and an O atom. Evolution of surface O abundance seems to be the key to 
the difference between the RE and MRE models. In the RE model the number of O atoms 
on a surface is determined by the balance between their accretion from the gas phase and 
consumption in O + HNO reaction. The same processes define O abundance in the MRE 
model as well, but only at t < 10^ years. After this time 0+HNO reaction slows down 
significantly due to modifying correction, and O abundance is controlled purely by accretion 
and desorption thereafter (as it is in the MC model at all times). As a result, OH and H2O 
abundances decrease, getting closer to the MC model prediction. A similar mechanism is at 
work in the case of ammonia as well. 



J^.2.2. Model H 

The only parameter region where a discrepancy occurs for Model H is located at n = 
2 ■ 10^ cm~^ and T = 25 K. Because of the low density, the degree of molecular complexity at 
these conditions is also low, and by 10^ years only 73 species have abundances higher than 
10-^2 at least in one of the calculations. Of these 73 species, 25 show disagreement between 
the MC calculations and the RE calculations. 

In the top row of Figure [8] we show the low density evolution of gas-phase water and 
molecular oxygen abundances. Both plots show a similar trend. In the MC calculations the 
abundance of the molecule stays nearly constant. At earlier times this steady-state behavior 
is reproduced in the RE run as well, but later the abundance grows and exceeds the MC 
abundance by more than an order of magnitude. It is interesting to note that the average 
O2 abundance in the MC run seems to decrease somewhat at later times, and this trend 
is reproduced by the RE model. The behavior of gas-phase abundances is a reflection of 
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grain-surface abundances. It must be kept in mind that surface water abundance is quite 
high in the H-RE model, reaching almost 10~^. Even though evaporation rate is not very 
high, it does increase gas-phase abundance up to a level of a few time 10~^°. Also, it is a 
low density model, so desorption is primarily caused by photons. 

In Figure [9] we show the evolution of H2O2 surface abundance and O2 gas-phase abun- 
dance in Model T-RE at n = 2-10^ cm~^ and the gas-phase ammonia abundance in Model H- 
RE at n = 2 ■ 10^ cm~^. These plots demonstrate that disagreement does not always mean 
simply higher or lower abundances. Different processes determine the H2O2 abundance at 
different times, and while the RE method is able to capture some of these processes, others 
are obviously missed. Our analysis shows that this behavior is related to the treatment 
of reaction with H atoms. Specifically, surface abundance of H2O2 is controlled by three 
reactions only. It is produced in a reaction H + O2H and destroyed by atomic hydrogen in 
reactions producing either H2O + OH or 0+2II + II2. (Note that in the considered model 
water ice is almost irreversible sink for O atoms as there are no surface reactions destroying 
water, and desorption is negligible.) The formation reaction is fast as it does not have an 
activation barrier. So, initially, as H and O atoms start to stick to a grain, surface abun- 
dance of H2O2 steadily grows. At the same time, there is not enough H atoms for either 
of destruction reactions. Then, in the MC model at around 1300 years the first destruction 
reaction (H + H2O2 — »■ H2O -|- OH), which has a slightly lower activation barrier than the 
other one (H + H2O2 — > O2H + H2), starts to "steal" some hydrogen atoms from the for- 
mation reaction, irreversibly removing one O atom per reaction from the H2O2 (re)formation 
process. This is the reason for the sharp fall-off of H2O2 abundance. On the other hand, 
in the RE model both destruction and formation reactions may occur simultaneously, which 
prevents the first destruction reaction from having such a dramatic effect. Note that the 
second reaction restores an O2H molecule which may react again with an H atom to restore 
an H2O2 molecule. 

H2O2 is not unique in this kind of behavior. There are some other surface species 
which show more or less similar trends, that is, almost perfect agreement at some times and 
quite noticeable disagreement at other times. These are some carbon chains (C2H2, C2N) 
or simpler molecules (CS, NO). Another equally dramatic example is represented by C2H2. 
The mechanism is similar, that is, it is related to the sequence of hydrogen additions which 
is treated differently in the MC model and the RE model. 

Evolution of the O2 abundance computed with the RE method follows the evolution 
computed with the MC model, but most of the time the RE curve is somewhat higher than 
the MC curve. Finally, the plot for ammonia shows the situation when the RE method 
correctly describes the average abundance evolution. However, the abundance predicted by 
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the MC method fluctuates so wildly, that at each particular t we can detect disagreement 
with high probability, using our formal order-of-magnitude criterion. 

5. Discussion 

In the previous sections we investigated the validity of the (modified) rate equation 
approach to grain surface chemistry under different physical conditions encompassing diffuse 
clouds, giant molecular clouds, infrared dark clouds etc. Evolution of the medium is studied 
with an extended gas-grain chemical model over a time period of 10^ years. We utilize the 
unified stochastic Monte Carlo approach, applied simultaneously to gas-phase and grain- 
surface reactions. The results are used to test the validity of conventional deterministic 
approaches. In general, differences in results obtained with deterministic and stochastic 
methods strongly depend on the adopted microphysical model of surface chemistry. In Model 
T, where tunneling for atomic and molecular hydrogen is permitted and diffusion/desorption 
energies ratio is low {E\^/Eiy = 0.3), discrepancy both for RE and MRE methods is very 
significant. At low temperatures (10 - 20 K) only abundances of surface species are discrepant 
by more than an order of magnitude. At moderate temperatures (25-40 K), due to active 
gas-grain interaction, incorrect treatment of grain surface chemistry becomes important for 
gas-phase abundances, too, and leads to dramatic decrease of overall agreement. At even 
higher temperatures agreement becomes better again due to limited surface chemistry. In 
Model H with no tunneling and high E\y/E-£, ratio, surface species are much less mobile, and 
the agreement between deterministic and stochastic methods is better. The only exception 
is low density (n(II) = 10^ cm^^) with moderate temperatures (25 - 30 K), where average 
residence time of species on grain surface is still lower than the average interval between 
reaction events. In general, the stochastic chemistry severely affects abundances and, thus, 
must be taken into account in chemical models of a warm and moderately dense medium. 

Modifications in rate equations, aimed at taking into account the accretion-limited 
regime, do not provide a significant improvement over 'canonical' rate equations in Model T. 
Close inspection of Figure [2] shows that at earlier times {t < 10^ years) modified rate equa- 
tions mostly produce results, which are rather less consistent with the 'exact' Monte Carlo 
solution, than results of the RE method. Only at later times {t ~ 10^ years) results of the 
MRE method become more consistent with the 'exact' solution than the RE results. In 
the previous section, we have already shown some examples of species, for which the MRE 
method actually worsens the agreement. 

Because at many of considered physical conditions results obtained with Monte Carlo 
method deviate significantly from those obtained with (modified) rate equations, the natural 
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question to ask is whether the inabihty of the RE method to treat these combinations of 
parameters might have caused mis-interpretation of observational data. To answer this 
question, one would need to consider an object, for which a large volume of observational 
data is available, so that simultaneous comparison for many molecules is possible. So far, 
there seems to be just one molecular cloud, which is studied in necessary details. This is a 
sub-region of the TMCl cloud where cyanopolyynes reach high concentrations, the so-called 
TMCl-CP. The physical cond itions in this object (T^ 10 K, nn ~ 10^ cm~'^) do fall into the 



range adopted in our model (IBenson fc Myerslll989l ). However, the TMCl-CP is probably 



too cold to show significant dependence of gas-phase abundances from surface processes. 

We may expect that the disagreement will be more significant in hotter objects. The 
range of physical conditions studied in this paper is quite narrow, as it has only been chosen 
to provide an initial view on the importance of stochastic treatment for surface chemistry. 
There are two directions which are worth to follow for further study. One of them is related 
to diffuse clouds. In this paper we most of the time assumed initial conditions in which all the 
hydrogen is molecular even at low density. Another simplification, which is, strictly speaking, 
not valid in low density medium, is equality of dust and gas temperature. So, it is interesting 
to check if stochastic effects work in predominantly atomic gas with different gas and dust 
temperatures. Another direction is chemistry at later stages of star formation, that is, in hot 
cores, hyper- and ultracompact HII regions, protoplanetary disks. In particular, in dark cold 
disk midplanes, which are poor in atomic hydrogen, stochastic effects may severely influence 
surface synthesis of hydrogenated species, like formaldehyde, which are later tran sferred by 



turbulent mixing in warmer disk areas and desorbed into gas phase (I Aikawal 120071 ). In these 
circumstances one would definitely need to use more accurate methods to model surface 
chemistry. Stochastic effects may be even more pronounced in the dynamically evolving 
medium, like during slow warm-up phase in hot cores or "corinos". Complex (organic) 
molecules will be among the most affected species. 

While the Monte Carlo method is the most direct approach, it is by no means practical. 
So, other methods are to be developed, in particular, the method of moment equations, if we 
want to understand deeper the chemical evolution in warm and dense astrophysical objects. 



6. Conclusions 

In this study, a gas-grain chemical model is presented, consisting of about 600 species 
and 6000 reactions with stochastic description of grain surface chemistry. For the first 
time both gas phase and grain surface reactions are simulated with a unified Monte Carlo 
approach. This unified model is used to test the validity of rate equations and modified rate 
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equations for a set of physical conditions, relevant to translucent clouds, dark dense cores, 
and infrared dark clouds. 

A comparison of results obtained with deterministic rate equations (RE) and modified 
rate equations (MRE) with results of Monte Carlo simulations was performed for this set 
of physical conditions and two different models of surface chemistry. We found that results 
obtained with both RE and MRE approaches sometimes significantly deviate from those 
obtained with Monte Carlo calculations in Model T, where surface species are highly mo- 
bile and most of the grain chemistry occurs in the accretion limited regime. While at low 
temperatures (10 K - 20 K) mainly RE and MRE abundances of surface species show signifi- 
cant discrepancies, at moderately high temperatures (25 K - 30 K) even abundant gas-phase 
species hke CO considerably deviate from the Monte Carlo results. At these temperatures 
grain surface chemistry is still very active and at the same time processes of accretion and 
desorption are very efficient. The gas-phase abundances of many species are influenced by 
surface chemistry that is correctly described in Monte Carlo simulations only. In Model 
H surface species are much less mobile than in Model T and the agreement between the 
(modified) rate equations approach and Monte Carlo simulations is almost perfect. The only 
parameter region where this discrepancy is significant is at a temperature of T = 25 K and 
at a low density n(H) = 2 • 10^ cm~^. On the whole, results of modified rate equations 
seem to be in somewhat closer agreement with Monte Carlo simulations than results of the 
conventional rate approach. But still, the agreement is far from being perfect, so modified 
rates should be used with care. We conclude that stochastic effects have a significant impact 
on chemical evolution of moderately warm medium and must be properly taken into account 
in corresponding theoretical models. In the studied parameter range, stochastic effects are, 
in general, not important at low (T ~ 10 K) and high (T ~ 50 K) temperatures. 
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Table 1: Adopted surface chemistry models 



Model 


Based on 


Tunneling 




T 


Hasee;awa et al. fl992) 


Yes 


0.3 


H 


Katz et al. fl999) 


No 


0.77 
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Fig. 1. — The Monte Carlo time step as a function of time in the model of the chemical 
evolution of a dense cloud. 
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Fig. 2. — Global agreement diagrams for Model T at three time moments. Percentage of 
species is shown in contour labels, for which abundances computed with the rate equations 
and with the Monte Carlo method differ by no more than an order of magnitude. Grayscale 
map with arbitrary contours is added for clarity. Darker color corresponds to worse agree- 
ment. Results are presented for conventional rate equations (top row) and for modified rate 
equations (bottom row). 
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Fig. 3. — Same as in Figure [2], but for Model H. 
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Fig. 4. — Evolution of H2CCCC gas-phase abundance at n = 2 • 10^ cm~^ and T — 10 K, 
simulated with the MC and RE methods (Model T). Dotted line indicates raw MC results, 
solid line shows MC results after averaging. Results of RE integration are shown with dashed 
line. 
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Fig. 5. — (Top row) Agreement diagrams for CO in Model T. Contour labels correspond to 
the logarithm of the ratio of gas-phase CO abundance in the MC run to that in the RE (top 
left) and MRE (top right) runs. (Bottom row) Gas-phase CO abundance evolution in the 
RE run (n = 2 ■ 10^ cm'^ and T = 30 K; bottom left) and in the MRE run (n = 2 ■ 10'' cm'^ 
and T = 25 K; bottom right). 
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Fig. 6. — (Top row) Agreement diagrams for methanol. Contour labels correspond to the 
logarithm of the ratio of gas-phase CH3OH abundance in the MC run to that in the RE (top 
left) and MRE (top right) runs. (Bottom row) Gas-phase methanol abundance evolution in 
the RE run (n = 2-10^ cm~^ and T = 25 K; bottom left) and in the MRE run (n = 2-10^ cm"^ 
and T = 25 K; bottom right). 
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Fig. 7. — Water and ammonia gas-phase abundance evolution in the RE run in Model T 
(n = 2 ■ 10^ cm-3 and T = 25 K; left) and in the MRE run (n = 2 ■ 10^ cm"^ and T = 25 K; 
right). 
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Fig. 8. — Water and molecular oxygen abundances at n = 2 ■ 10^ cm ^ and T = 25K in 
Model H. 
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Fig. 9. — Abundance evolution for selected species. Models and physical conditions are 
indicated on the diagrams. 



